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ABSTRACT 

The process of exploring the behavior of nonlinear, dynamical systems can be a time- 
consuming and tedious process. In this thesis, I have written a program which automates 
much of the work of an experimental dynamicist. In particular, the program automat¬ 
ically characterizes the behavior of any driven, nonlinear, electrical circuit exhibiting 
interesting behavior below the 10 Mhz range. In order to accomplish this task, the pro¬ 
gram can autonomously select interesting input parameters, drive the circuit, measure 
its response, perform a set of numeric computations on the measured data, interpret 
the results and decompose the circuit’s parameter space into regions of qualitatively dis¬ 
tinct behavior. The output is a two-dimensional portrait summarizing the high-level, 
qualitative behavior of the nonlinear circuit for every point in the graph as well as an 
accompanying textual explanation describing any interesting patterns observed in the 
diagram. In addition to the graph and the text, the program generates a symbolic de¬ 
scription of the circuit’s behavior. This intermediate data structure can then be passed 
onto other programs for further analysis. 


This report is a revised version of a thesis submitted on May 11, 1989 to the Depart¬ 
ment of Electrical Engineering and Computer Science at the Massachusetts Institute of 
Technology in partial fulfillment of the requirements for the degree of Master of Science. 
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Chapter 1 


Introduction 


1.1 Goals 

This thesis is motivated by the desire to automatically measure the behavior of complex, 
dynamical systems and to generate high-level, qualitative interpretations of their behav¬ 
ior. In this thesis, I have put together a system capable of automatically measuring the 
behavior of electrical circuits and performing this high-level interpretation of its behav¬ 
ior. The system takes advantage of information provided by sophisticated measurement 
instruments and numerical software packages to aid in the identification and analysis of 
the circuit’s behavior. 


1.2 Background and Motivation of Thesis 

Physicists, mathematicians, and engineers have spent a great deal of time studying the 
behavior of nonlinear, dynamical systems. They have written numerous papers describ¬ 
ing each system’s behavior, focusing in particular on the transition route from stability 
to chaos[DBHL82, Hub83, CYY86, Lin81, KC86, TPJ82, UA81, KKC, MCK85, Sha81]. 
These papers typically describe interesting behaviors or patterns exhibited by the dynam¬ 
ical system and the parameter values at which they occur. For example, in Evidence for 
Universal Chaotic Behavior of a Driven Nonlinear Oscillator[ TPJ82], Testa, Pdrez, and 
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Jefferies discuss the parameter values at which they observe the period-doubling 1 route 
to chaos in their driven, nonlinear oscillator. Similarly in The Devil’s 5fairawe[KKC], 
Kennedy et. al. trace the progression of the onset of chaos 2 in their transistor circuit; 
and Chua, in the Devil’s Staircase Route to Chaos in a Nonlinear C'treu»f[CYY86], ex¬ 
plores the parameter values at which his oscillator circuit exhibits subharmonics of all 
orders from one to infinity. In fact, the first paragraph describing Chua’s experimental 
results in the Devil’s Staircase paper reads, 

As we decrease the frequency /, from 20 khz to 500 khz, with all other 
parameters held fixed, we observe subharmonics of all orders from 2 to 33. 

Each subharmonic is found to persist over a limited range of input frequency 
/». Each [point] in the [parameter-space] graph can be interpreted as a syn¬ 
chronization state between the input frequency and some submultiple of the 
circuit’s natural frequency, i.e. the oscillation frequency when the input sig¬ 
nal is set to zero. The [area of each region] can therefore be interpreted as 
the “locking range.” [The greater the area; the greater the locking range.] 

Aside from providing textual descriptions, investigators have also generated parameter- 
space graphs which summarize the behavior of a complex, nonlinear system for each 
point in the two-dimensional grid. For example, Linsay, in his paper Quasiperiodicity 
and Chaos in a System with Three Competing Jf'Veguenctes[Lin.88], presents a series of. 
parameter-space diagrams one of which is shown in figure 1.1, and similarly Hayashi 
makes use of parameter-space diagrams as shown in figure 1.2[Hay64]. 

Exploring the behavior of a dynamical system over the range of input parameters 
clearly provides useful information. Mechanical engineers, for instance, spend much of 
their time studying the behavior of the mechanical systems they have designed. They 

1 A nonlinear system is said to exhibit the period-doubling route to chaos if its behavior progresses 
from subharmonics of order l-»2-»4-»8-» 16-» 32chaos. 

3 Chaos is a commonly observed phenomenon in a wide range of physical systems such as nonlinear 
electrical circuits, control systems with nonlinear restoring forces, thermal convection in fluids, vibrations 
in the panels of a supersonic aircraft or rocket, the swings of a pendulum with a vibrating pivot, the 
rotations of Hyperion, a moon of Jupiter, buckled elastic beams, certain chemical reactions, propagation 
of light through nonlinear optical devices, the ventricular fibrillation of a heart beat, the firing of neurons 
in marine mollusk and even the dripping of a leaky kitchen iaucet[Moo87]. 
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Figure 1.1: Example of a parameter-space graph from Unsay. Points at which the circuit 
exhibits a periodic response are shown by dots. The white regions indicate quasiperiodic 
behavior, and the black squares represent chaos. 

are typically concerned with locating regions of stability and instability in their devices. 
Having parameter-space diagrams such as the ones shown in figure 1.1 and figure 1.2 
allows the investigator to locate the regions of interest at a glance. 

Referring to figure 1.2, we find that at input values of / = 0.5 and A = 10 (the point 
labelled 1), the nonlinear system exhibits harmonic behavior whereas point 2, at / = 2 
and A = 0.25, lies on the ^-harmonic and the ^-harmonic border. The diagram also 
decomposes the system’s parameter space into regions of qualitatively distinct behavior. 
For instance, for values of A and / near the origin, the dynamical system displays brief 
occurrences of fifth-order and third-order harmonic behavior. For values of A and / far 
away from the origin, the system exhibits higher order subharmonic behavior such as 
fifth and sixth. 

Exploring the behavior of a complex, nonlinear system and producing a comprehen¬ 
sive parameter-space diagram such as the ones shown in figure 1.1 and figure 1.2 can be a 
time consuming and tedious process[AS87]. There exist many techniques to help identify 
a system’s behavior, including power-spectrum analysis, tests for periodicity, and the 
stroboscopic technique. However, the investigator must still repeatedly select interesting 
values for parameters and initial conditions, run the test, interpret the results, classify 
the behavior and record it. This must be done for every point in a two-dimensional 
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Figure 1.2: Example of a parameter-space graph from Hayashi. 


parameter-space. Moreover, in cases where the behavior is not easily identifiable, the 
investigator must run additional tests to gain different perspectives of the system’s be¬ 
havior. Using the computer to automate this process can greatly aid the investigator in 
his initial exploration of the system’s behavior. 

In this thesis I describe a program that can automatically produce a two-dimensional 
parameter-space diagram characterizing the behavior of a dynamical system. It can 
set up, run, and interpret measurements for each point in the parameter space. The 
program can identify four types of behavior: harmonics, subharmonics, quasiperiodics, 
and chaos. 3 In addition to the graph, the program also produces a high-level, textual 
explanation to describe any interesting behaviors or patterns in the parameter space 
of the system. 4 In particular, the textual explanation draws attention to large and 
small regions of uniform behavior, large regions of chaotic behavior, frequency-locking 


3 Precise definitions of these terms and examples will be given in subsequent chapters. 

4 In the past, scientists like Ueda and Linsay have written programs to automate parameter-space 
graph generation[UA81, Lin88]. However, their programs did not interpret the high-level behavior in 
the graph. Instead, they left it up to the user of the program to notice any interesting patterns such as 
period-doubling cascades. In fact, many of the published papers in the literature of nonlinear dynamics 
spend a considerable amount of time discussing the investigator’s high-level interpretations of the circuit’s 
behavior. Chna’s quote on page 10 provides a brief example. The program in this thesis goes beyond the 
work of Linsay et. al. by automatically describing the interesting, high-level, qualitative characteristics 
of the parameter-space portrait. With such a capability, the program in this thesis could potentially 
write parts of many of the published papers in the literature of nonlinear dynamics. 
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ranges, 5 period-doubling cascades and period-adding cascades. 6 It also explains how the 
behavior of the system varies as the user varies a single input parameter, and it describes 
the bifurcation curves 7 which border each region of uniform behavior. This is the kind 
of qualitative, high-level information an investigator seeks when studying the behavior 
of complex, nonlinear systems experimentally. 

Finally, after computing the graph and generating the accompanying textual expla¬ 
nation, the program also returns a symbolic description of the circuit’s behavior in the 
form shown in figure 1.3. This symbolic descriptor, containing information about regions 
of uniform behavior and boundaries between regions, might then be used by another pro¬ 
gram to further analyze the behavior of the circuit. For instance, a program might take 
the symbolic description as input, and it might suggest to the user all the regions in the 
parameter space in which the device operates in a stable manner. Or it might search 
the graph for other interesting patterns such as saddle-node, pitchfork, flip or Hopf 
bifurcations[TS86]. Or it might locate all regions of instability and transition and warn 
the user against operating near these regions. The program is capable of offering such 
advice because information about the circuit’s behavior is summarized in the symbolic 
description. 

1.3 Sample Run of the Analyzer Program 

Figure 1.4 provides a sample parameter-space graph generated by running the program 
on the forced negative-resistance oscillator circuit[UA81]. Also, accompanying the graph 
is a text explanation describing the high-level, qualitative properties of the behavior of 
the circuit. To start the program running, the user needs to type the following input 

5 Frequency-locking ranges are regions in a nonlinear system’s parameter-space where the behavior 
remains constant despite variations to an input parameter. For instance, if we increase the driving 
frequency of a given circuit from 30000 hertz to 60000 hertz and its response remains fixed at subharmonic 
order 1, then the behavior of the circuit is said to be “locked” in the range from 30000 hertz to 60000 
hertz. 

* A nonlinear system is said to exhibit a period-adding cascade if its behavior goes from subharmonics 
of order 1—►2—>3—>4—>5—►6—.7—.8—chaos[CYY86]. Compare this definition with the 
definition of period doubling as defined in an earlier footnote on page 1. 

7 A bifurcation is a qualitative change in the organization of an equilibrium state as a system parameter 
is varied. 
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( ;; lumped regions of uniform behavior 
#(((SUBHARMONIC 3) ;; first region of uniform behavior 

(((VOLTAGE 9) (DRIVE-FREQUENCY 47000)) ;; points in region 

((VOLTAGE 8.5) (DRIVE-FREQUENCY 47000)))) 

((CHAOTIC) ;; second region of uniform behavior 

(((VOLTAGE 9) (DRIVE-FREQUENCY 47500)) ;; points in region 

((VOLTAGE 8.5) (DRIVE-FREQUENCY 47500)))) ) 


;; boundary information 
(((SUBHARMONIC 3 "Rl") 
((CHAOTIC "R2"))) 

((CHAOTIC "R2") 
((SUBHARMONIC 3 "Rl")))) ) 


Figure 1.3: Sample symbolic description returned by the program. The first part of 
the symbolic descriptor provides information concerning the lumped regions of uniform 
behavior. The second part presents boundary information. 

command: 

General Format: 

(run-point-by-point-2D-analysis circuit-name test-number 

starting-volt ending-volt volt-units 
starting-freq ending-freq freq-units 
volt-delta freq-delta) 

Specific Example: 

(run-point-by-point-2D-analysis ’neg-resis 4 

1.5 6.5 v 
40e3 70e3 hz 
.25 le3) 

where the first argument specifies the name of the circuit under investigation. The second 
indicates the test number. The third, fourth, and fifth arguments specify the starting 
voltage, ending voltage and voltage units respectively. Similarly, the sixth, seventh, and 
eighth arguments pertain to the starting frequency, ending frequency and frequency units 
respectively. The last two parameter values indicate the incremental step size for voltage 
and frequency. Thus, in this particular example, the program will analyze the behavior 
of the forced negative-resistance oscillator over the frequency range of 40000 hertz to 
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70000 hertz at increments of 1000 hertz and over the voltage range of 1.5 volts to 6.5 
volts at increments of 0.25 volts. 

Given this information, the program automatically drives the oscillator circuit with 
a sinusoid of the appropriate amplitude and frequency, sends the appropriate commands 
to a waveform recorder to record the signal, classifies the response, and produces the 
parameter-space diagram shown in figure 1.4 along with its accompanying high-level tex¬ 
tual explanation shown in figure 1.6. Notice how the textual explanation draws attention 
to the large and small regions of uniform behavior, the transitions between regions, and 
any special patterns such as period-doubling cascades. Notice also how the program, 
like a human, correctly identifies three main regions of uniform behavior—first-order 
subharmonic, second-order subharmonic, and chaos— in a parameter space that really 
consists of nine or more regions. This illustrates the program’s ability to abstract the 
regions of uniform behavior, even in a messy graph. In other words, the program is able 
to take a global enough perspective so as to not let the small, isolated, and sporadic 
occurrences of third, fourth, fifth, seventh and ninth order subharmonics distract from 
the dominant first-order subharmonic, second-order subharmonic and chaotic behaviors. 
As will be discussed in chapter 5 of this thesis, the program accomplishes this abstraction 
by cleaning up or transforming the “messy” parameter space in figure 1.4 to the new 
graph shown in figure 1.5. 

Thus, the parameter-space graph, along with the textual description, gives the inves¬ 
tigator a quick and thorough understanding of the system’s behavior, drawing attention 
to the regions of qualitatively interesting behavior. After r unning the program once for 
a first-pass overview, the program can be easily run again, over a much smaller region, 
to obtain a more detailed picture. Each successive run allows the investigator to zoom 
in on the region of interest. 

1.4 Implications of Thesis Work 

The integration of sophisticated measurement instruments with powerful numerical pack¬ 
ages and software programs, which interpret the results of physical measurements and 
numeric computations, suggests new possibilities for test and measurement. 
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Figure 1.4: Sample parameter-space graph from the forced negative-resistance oscillator 
circuit. The A along the y-axis indicates the amplitude of the driving sinusoid, while the 
/ along the z-axis specifies the frequency of the driving sinusoid. 
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Figure 1.5: Cleaned up parameter space without blemishes. 
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HIGH-LEVEL INTERPRETATION OF PARAMETER-SPACE DIAGRAM 

The parameter-space diagram characterizes the behavior of the forced 
negative-resistance oscillator as the frequency and amplitude of the 
driving sinusoid are varied from 40000 to 70000 hertz and from 1.5 to 6.5 
volts respectively. For amplitudes between 1.5 and 6.5 volts and 
frequencies between 40000 to 56000 hertz, the system exhibits 
predominantly first order subharmonic behavior, while for amplitudes 
between 1.5 and 5.75 volts and frequencies between 51000 to 68000 hertz, 
the system exhibits predominantly chaotic behavior, while for amplitudes 
between 1.5 and 6.5 volts and frequencies between 55000 to 70000 hertz, 
the system exhibits predominantly second order subharmonic behavior. 

Aside from these luge areas of uniform behavior, the system passes briefly 
through subharmonics of order 7, 5, 9, and 4 at isolated points throughout 
the parameter-space graph. 

The large region of first order subharmonic behavior is bounded by the 
bifurcation curve(s) 1-C, and 1-2, where C is a region consisting of 
primarily chaotic behavior, but littered with sporadic occurrences of 
other kinds of behavior. Similarly, the luge region of chaotic behavior 
is bounded by the bifurcation curve(s) C-2, and C-l, while the region of 
second order subharmonic behavior is bounded by the bifurcation curve(s) 
2-C, and 2-1. 

Finally, based upon the data gathered so fu, the forced 
negative-resistance oscillator shows no cleu signs of period-doubling 
bifurcations as described by Feigenbaum and others. 

Figure 1.6: Text generated by the high-level interpreter. In the second paragraph, the 
notation 1-C indicates a boundary with first-order subharmonic behavior on one side 
and chaotic behavior on the other. Similarly, the notation 1-2 indicates a boundary with 
first-order subharmonic behavior on one side and second-order subharmonic behavior on 
the other side. 
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Currently instruments play a rather isolated and limited role in the design and ex¬ 
ploration stages. Designers or investigators working in most programming environments 
do not usually use or have access to a wide range of powerful measurement instruments. 
Even if they do have access to some instruments from their programming environments, 
the control that they have over an instrument and the results returned by an instrument 
tend to be fairly low-level. For example, a user might be able to send amplitude and 
frequency values to set up a function generator, or perhaps a dc-voltage level to a volt 
meter. Or a dynamic signal analyzer might return several thousand data points as a 
result of computing the fast Fourier transform (FFT) of a signal. But, interpreting this 
stream of data points to arrive at some sort of high-level understanding of the behavior 
still requires a fair amount of human effort. 

Aside from the issue of instruments, scientists have typically used computers almost 
solely for purposes of numeric computations such as convolutions and autocorrelations. 
The results of such computations are usually in the form of large bodies of numerical 
data. Although the computer is very good at performing numeric computations, it is 
not very good at interpreting its qualitative content. It is currently largely up to the 
investigator to scan through reams of data searching for relevant information. 

In this thesis, I have taken advantage of both sophisticated measurement instruments 
and powerful numerical packages. Yet, on top of these two tools, I have written a program 
which interprets the results of the numeric computations and the physical measurements. 
As a result, we now have a system with all the power of a whole range of sophisticated 
measurement instruments, a whole set of numerical packages, and a high-level interpreter, 
capable of producing meaningful and informative, qualitative, descriptions of complex, 
nonlinear systems. No longer are scientists limited to low-level, control-the-volt-meter 
type commands. They now have access to more powerful, more intelligent tools, capable 
of generating such high-level statements as, “The forced negative-resistance oscillator 
undergoes a period-doubling cascade as the driving frequency is increased from 26 khz 
to 78 khz while the amplitude is fixed at 5.25 volts.” 
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1.5 Organization of Thesis 

Chapter 2, presents a brief overview of the components within this thesis and how they 
interrelate. Chapters 3 and 4 discuss the details of how the analysis progr am classifies the 
behavior of the recorded waveform. In particular, chapter 3 focuses on the use of power- 
spectrum data for behavior identification, while chapter 4 focuses on the use of time-signal 
information. In chapter 5, I describe the high-level interpreter which is responsible for 
recognizing interesting patterns in the parameter-space graph and 
level, textual explanation which accompanies the graph. Finally, 
with an analysis of the contributions of this thesis and suggestions 


generating the high- 
chapter 6 concludes 
for future work. 
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Chapter 2 


Overview 


This chapter presents the experimental set up of the system as well as a brief description 
of each component within this system. As illustrated in figure 2.1, the system consists 
of three main modules: the Observation Equipment, the Measurement Instruments, and 
the Software Module. 

2.1 The Nonlinear Circuit 

Before discussing these three components, I would first like to describe the types of 
circuits the program can handle. Any driven, nonlinear, electrical circuit exhibiting 
interesting behavior below the 10 Mhz range can be connected to this system for analysis. 1 
Moreover, the current system can only have two free parameters: driving amplitude and 
frequency. The nonlinear circuit shown in figure 2.1 is Ueda’s forced negative-resistance 
oscillator circuit[UA81]. I specifically chose to study this circuit because it provides 
a rich source of interesting, nonlinear behavior- subharmonics of orders 1 through 19, 
quasiperiodics, and even chaos. This oscillator is driven by a sinusoid generated by an 
HP 3325a function generator, and its response is measured over the capacitor and input 
into the HP 5182a waveform recorder and the oscilloscope. Its nonlinear behavior stems 
from the nonlinear resistor parallel to the capacitor. Figure 2.2 illustrates the schematic 

The 10 Mhz constraint is a limitation set by the HP 5182a Waveform Recorder/Generator. Refer to 
the specification sheet in the Operating Manual for the HP 5182a.[51884]. 
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Figure 2.1: Overall layout of the system. The measurement instrument component drives 
the nonlinear circuit with a sinusoid and records the circuit’s response. The software 
module analyzes and interprets the results of the measurements, and the observation 
equipment allows the user to visually monitor the behavior of the circuit. 
















Figure 2.2: Nonlinear resistor in the forced negative-resistance oscillator. 

of the nonlinear resistor in the forced negative-resistance oscillator. 

In addition to the forced negative-resistance oscillator, we have constructed sev¬ 
eral other circuits exhibiting interesting nonlinear behavior— Chua’s devil’s staircase 
circuit[KKC], Matsumoto’s double scroll[MCK85], and a simple, series RLC circuit where 
the nonlinear element is a varactor diode whose capacitance is a function of voltage. 
These circuits were connected to the system and used to test the robustness of the anal¬ 
ysis program, which was developed primarily using Ueda’s forced negative-resistance 
oscillator circuit. 


2.2 The Measurement Instruments 

As mentioned earlier, the system uses a set of sophisticated measurement instruments. 
In particular, an HP 5182a waveform recorder /generator is used to record the response 
waveform of the circuit under investigation and an HP 3325a synthesizer function gener¬ 
ator to produce the sine wave that drives the circuit. For circuits exhibiting interesting 
behavior under the 100 khz region, the HP 3562a dynamic signal analyzer can also be 
hooked up to the system for measurement and analysis. However, since most of the in¬ 
teresting behaviors occur at the higher frequencies, the dynamic signal analyzer has been 
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of limited use to us. All instruments axe controlled remotely by an HP 320 computer via 
HP-IB 2 commands. 

2.3 The Observation Instrument 

The observation component of the system does nothing to help the analysis program clas¬ 
sify the nonlinear circuit’s behavior. Its only purpose is to let the user visually monitor 
in real time the signal and power spectrum of the response. As shown in figure 2.1, the 
driving sinusoid of the circuit is fed into channel A of the oscilloscope while its response 
is fed into channel B of the scope. In addition, the output of an HP 8007a pulse generator 
is connected to the 2 -input channel at the back of the scope. When the scope is put into 
x — y mode, we get a Lissajous figure. By counting the number of crests in the Lissajous 
figure, it is possible to determine the order of the subhaxmonic[CYY86]. This technique, 
however, is only practical if the subharmonic order is small enough so the eye can dis¬ 
tinguish one crest from another on the oscilloscope’s screen. An easier way to determine 
the subharmonic order by inspection is to strobe the time signal with the pulse generator 
of the appropriate frequency. If a response is first-order subharmonic, it will display a 
single highlighted dot in the Lissajous plot. If it is second order, it will have two dots 
in the Lissajous plot. If the response is chaotic or some very high-order subharmonic 
response, the Lissajous plot will show a blurring of many, perhaps an infinite number, of 
dots. Counting the dots is an easy way to determine order subharmonic. 

In addition to the oscilloscope, I have input the circuit’s response into an HP 3562a 
dynamic signal analyzer. The 3562a has the ability to compute and display on its screen 
the power-spectrum response of a signal in real time. Both the oscilloscope and the 
dynamic signal analyzer provide useful tools to monitor the behavior of the response. 

2.4 The Software Module 

As shown in figure 2.1, the software module consists of the analyzer program, the 6.003 
signal-processing package, and the communication module. The code for all three com- 

3 Hewlett-Packard Interface Bus 
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ponents resides on the HP 320 computer. 


2.4.1 The Analyzer Program 

The analyzer program is the heart of the system. It is responsible for setting up the test 
runs, determining the appropriate parameter values for proper or optimal functioning of 
the instruments, and sending the appropriate HP-IB commands to the function generator 
to drive the circuit and to the waveform recorder to record the circuit’s response. Once 
the response has been recorded, the analyzer program transfers the recorded waveform 
from the waveform recorder to the computer, computes the power spectrum for each 
response, and classifies the behavior of the circuit into one of four categories: harmonic, 
subharmonic, quasiperiodic, or chaotic. After having identified the behavior of the circuit 
for every point in the two-dimensional grid, the analyzer program describes the high-level, 
qualitative characteristics of the nonlinear system. 

Figure 2.3 illustrates the block diagram of the analyzer program. The inputs to the 
program are data points representing the signal captured from the nonlinear circuit by the 
waveform recorder, and the outputs include a parameter-space graph, a high-level textual 
explanation describing the circuit’s qualitative characteristics, and a symbolic description 
to be passed to other programs. The behavior classification module of the system, shown 
in the top portion of figure 2.3, identifies a system’s behavior by examining both the 
power-spectrum data and the time-signal data. The high-level interpreter component 
of the system consists of all the boxes below the resolver in figure 2.3. Its purpose is 
to generate the qualitative text explanation. I will discuss the behavior classification 
portion of the analyzer program first. This consists of two parts: The power-spectrum 
analyzer and the time-signal analyzer. 

The Power-spectrum analyzer 

Once the recorded waveform has been down-loaded onto the computer, the power- 
spectrum analyzer computes its power spectrum. This power spectrum provides crucial 
information for the classification of the circuit’s behavior. Based upon the location of 
the peaks in the spectrum, the program can categorize the behavior of the circuit into 
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one of four categories: harmonic, subharmonic, quasiperiodic, or chaotic. Chapter 3 
provides an in-depth discussion of the peak detection and the subsequent classification 
of behavior. 

The technique of identifying a system’s behavior by searching for peaks in the power 
spectrum proves effective when the noise level in the spectrum is low. However, as 
the noise level in the spectrum increases, peak detection becomes more difficult and the 
accuracy of the power-spectrum analyzer’s conclusion degrades. To remedy this problem, 
the program looks at additional information contained in the time signal to confirm or 
disconfirm the analysis of the power-spectrum analyzer. By examining the behavior 
from two perspectives, the program increases the likelihood of at arriving at the correct 
conclusion. 

The Time-signal analyzer 

As mentioned in the previous section, the time signal also plays an important role in 
the identification of the system’s behavior. Referring to figure 2.3, we see that the time- 
signal analyzer strobes the response waveform at multiples of the period of the driving 
sinusoid. If the response repeats itself with every strobe, it exhibits harmonic behavior. 
If the response repeats itself every other time, it exhibits second-order subharmonic 
behavior; if every third time, third-order subharmonic response, etc. As the number of 
strobes needed for the signal to repeat itself approaches infinity, so too does the order 
subharmonic. Subharmonic responses whose order approaches infinity axe essentially 
chaotic since by definition a chaotic signal does not ever repeat itself. Chapter 4 discusses 
in greater depth the details of the time-si gn a.1 analyzer. 

The Resolver 

The conclusion produced by the power-spectrum analyzer and the conclusion of the time- 
signal analyzer should concur. However, if they do not, the resolver module, as shown 
in figure 2.3, has the responsibility of resolving the conflict. To decide which one of the 
two conclusions to accept, the resolver looks at contextual information. If, for example, 
the power-spectrum analyzer identified the behavior of a particular point as CHAOTIC, 
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and the time-signal analyzer classified the behavior of the same point as SUBHARMONIC 
ORDER 2, the resolver will look at the neighboring points to the east and to the south. 
If one or both are SUBHARMONIC ORDER 2, then the resolver will assume that, the power 
spectrum was extraordinarily noisy and will reject the CHAOTIC analysis in favor of the 
SUBHARMONIC ORDER 2 analysis. Chapter 4 presents the resolver in greater detail. 

The High-level interpreter 

Once the program has classified the behavior of the circuit for every combination of 
input parameters, it produces a two-dimensional parameter-space graph and passes the 
result to the high-level interpreter. The high-level interpreter scans the parameter-space 
diagram, first cleaning up any small, isolated occurrences of sporadic behavior in a messy 
graph, then searching for the boundaries between regions of uniform behavior. Next, 
the high-level interpreter generates the textual explanations describing any interesting 
patterns observed in the graph. Knowledge about these patterns is contained in an easily 
extensible knowledge base. Currently the high-level interpreter can recognize only period- 
doubling cascades, period-adding cascades and frequency-locking ranges. However, as 
we add more patterns to the knowledge base, the high-level interpreter will be able to 
recognize more sophisticated patterns. In Chapter 5, I present some other commonly- 
observed, well-studied patterns in nonlinear dynamics which we may want to include. 

2.4.2 The 6.003 Signal-processing Package 

Although the analyzer program makes up a large portion of the software module, the 
6.003 signal-processing software also has an important function. This signal-processing 
environment, developed at M.I.T. for the Circuits, Signals, and Systems course, provides 
a variety of tools for signal analysis, such as fast Fourier transforms, auto-correlations, 
chirp ^-transforms, and convolutions. The primary purpose for including this signal¬ 
processing package in the system is to take advantage of its FFT and its autocorrelation 
procedures. 
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2.4.3 The Communication Module 


The communication module is responsible for translating a high-level Scheme command 
into the appropriate sequence of HP-IB commands specific to the particular instrument. 
For example, the analyzer program contains such high-level commands as: 

(get-time-signal ... ) 

By instrument standards, this is a fairly high-level command. In order to successfully 
record a signal, the HP 5182a waveform recorder needs a plethora of additional informa¬ 
tion about the range of the signal, the channel of the incoming signal, whether to AC 
couple the signal, whether to select automatic sweep mode for triggering, the value for 
the time base, the memory record length, etc. To make matters worse, each instrument 
has a different set of HP-IB commands which the user must learn if he wishes to commu¬ 
nicate remotely with the instrument. Ideally the user of the system should not have to 
worry about such details every time he wishes to perform a simple task such as record¬ 
ing a signal. The communication module is smart enough to translate such high-level 
commands as get-time-signal into a sequence of HP-IB commands, with optimally 
chosen parameter values, to perform the appropriate task. 3 For example, the following 
code shows the commands generated by the communication module to configure the 
waveform recorder for get-time-signal: 

>> Select Recording mode as opposed to Generating Mode for the HP 
;; 5182a Waveform Recorder/Generator. 

(write-string-to-analyzer "RC" 5182a-device-file) 

;; Auto set to certain default values. See Page 3-52 of 5182a 
;; Operating Manual. 

(write-string-to-analyzer "AU" 5182a-device-file) 

;; Select channel A as input. 

(write-string-to-analyzer "CHI" 5182a-device-file) 

;; Select Range for channel A. 

3 By optimally chosen parameter values I mean the program knows how to, for example, choose the 
sampling frequency so as to avoid aliasing. More specifically, it chooses the sampling rate so that it is 
sufficiently below the Nyquist rate. 
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(write-string-to-analyzer "AR" 5182a-device-file) 
(write-string-to-analyzer (number->string 500e-3) 

5182a-device-file) 

;; Select AC coupling for channel A. 

(write-strings-to-analyzer 5182a-device-file "AC1") 

;; Select Sweeping Arming for Trigger to be Automatic Sweep Node, 
(write-strings-to-analyzer 5182a-device-file "SA2") 

>* Set up value for time base. Note units of time-per-sample 
;; is in seconds. 

(write-string-to-analyzer "MM" 5182a-device-file) 
(write-string-to-analyzer (number->string time-per-sample) 

5182a-device-file) 

;; Set up memory record length. 

(write-string-to-analyzer "LE" 5182a-device-file) 
(write-string-to-analyzer (number->string memory-record-length) 

5i82a-device-file) 

(write-string-to-analyzer "L01" 5182a-device-file)) 
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Chapter 3 


The Power-spectrum Analyzer 


The purpose of the power-spectrum analyzer module is to classify the system’s behavior 
as chaotic, quasiperiodic, subharmonic or harmonic, based upon information presented 
in the power spectrum. It achieves this by breaking the task into three main parts: clean 
up spectrum, detect peaks, and classify behavior. 

3.1 Cleaning up the power spectrum 

All real power spectra, including those of perfectly harmonic and subharmonic responses, 
have a certain amount of low-level noise scattered throughout the spectrum. Therefore, 
before passing the spectrum through the peak detector, the power-spectrum analyzer 
tries to identify and then to minimize the effects of noise. In figure 3.1(b), I have plotted 
the power spectrum of a periodic response on a logarithmic scale in order to magnify the 
presence of the low-level noise in the spectrum. 

This noise may be attributed to several factors including: 

• Quantization Error due to the analog to digital conversion of the signal by the 
waveform recorder, 

• Arithmetic Error in the FFT algorithm, and 

• Spectral Leakage due to a finite-length representation of an infinite-length input 
signal. 

Quantization error results from the representation of continuous signal amplitudes 
by a fixed number of digital levels within a predefined maximum and minimum analog 
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a) 


b) 



Figure 3.1: Power-spectrum plot to magnify low-level noise: (a) on a linear scale; (b) on 
a logarithmic scale. 

range. Each value of the continuous signal, x(n), must be encoded using B bits to 
obtain the quantized output, x^(n). The quantization error stems from rounding the 
input value, x(n), to the nearest quantization level, x^(n). Figure 3.2 illustrates how 
the input x(n) = 0.450 gets rounded to the nearest quanta, x^(n) = 0.500. 

Thus, the error due to quantization can be represented as 

e(n) = |x(n) - x 9 (n)| (3.1) 

It is this e(n) term which contributes to the low-level noise observed in the power spec¬ 
trum. 

In general, the quantization error is a function of B, the number of bits in the A/D 
converter and the distribution of the input signal over the range of the A/D converter. 
The greater the number of bits, the smaller each quantum and the less the error. 

The peak-detection module minimizes the quantization error by appropriately setting 
certain parameters on the HP 5182a waveform recorder. The HP5182a with its 10 bits 
of accuracy, 1 ranges from -511 to +512, a span of 2 10 = 1024 = 512 - (-512) points. 

1 For more detailed specification information on the HP5182a waveform recorder, refer to the Operating 
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Figure 3.2: Encoding Process Related to A/D Conversion with B = 4 and q = 0.125 
[AN83] 


The program automatically scales the amplitude of the input signal so that the incoming 
waveform makes optimal use of the HP5182a’s full range. In other words, regardless of the 
amplitude of the incoming signal, the program scales the signal so that it occupies nearly 
the full range of the HP5182a from -512 to +512, thereby minimizing the quantization 
error. 

Spectral leakage results when Fourier transform algorithms try to approximate in¬ 
finitely long signals with finite-length signals. In practice, one typically has access to only 
a finite portion of an infinite signal. Computing the fast Fourier transform (FFT) of a 
finite-length signal has the same eflect as taking a finite piece of the signal, replicating 
it and computing its Fourier transform. If, however, the size of the finite segment does 
not contain an integral number of cycles, a discontinuity in the signal results, as shown 
in Figure 3.3[DLH88]. 

The discontinuity in the signal gives rise to the leakage of sidelobes (or spectral con- 
and Programming Manual[51884]. 
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Figure 3.3: Discontinuity resulting from periodic extension of signal by DFT process: 
(a) continuous signal; (b) signal showing discontinuity due to DFT. 

tributions) across the entire frequency set. (See Figure 3.4.) These high sidelobe levels 
in the spectrum can result in false peak detections. A common way of reducing sidelobe 
levels is by multiplying the signal by a weighting function[AN83, DLH88, OS75, Har78]. 
Weighting functions reduce the contribution of the samples near the endpoints, and there¬ 
fore reduce the discontinuity and its effect on the frequency response. Extensive evalu¬ 
ations comparing the advantages and disadvantages of the various windowing functions 
ranging from Kaiser-Bessel windows to Gaussian windows have been performed[Har78]. 
I have chosen to multiply the signal by a Hanning window of order one before taking the 
FFT. Figure 3.5 illustrates the the power spectra of a windowed input signal. Compare 
the differences between figure 3.4 and figure 3.5, and notice the reduction in height of 
the sidelobes due to spectral leakage. 

3.2 Peak-detection Component 

As mentioned in chapter 1, the key to classifying the behavior of a reponse using power- 
spectrum data lies in identifying the peaks and reco gnizing the relationships among their 
locations. Without an accurate peak detector, the power-spectrum analyzer module has 
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Figure 3.4: Power spectrum of a non-windowed signal, (a) The top graph illustrates 
the full spectrum on a Unear scale, (b) The bottom graph shows a magnified version 
of (a) around the peak to emphasize the effect of the sidelobes on the low-level noise 
throughout the spectrum. 



Figure 3.5: Power spectrum of a windowed signal, (a) The top graph illustrates the 
full spectrum on a Unear scale, (b) The bottom graph shows a magnified version of (a) 
around the peak. 
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very little chance of accurately classifying the behavior of a system. The purpose of the 
program’s peak-detection component is to sequentially check each point in the power 
spectrum to see which points possess the traits of a peak. Once it recognizes a peak, this 
module records the peak’s height and location so the classification module can classify 
the recorded waveform’s behavior. 

3.2.1 Characteristics of a Peak 

Each point in the power spectrum must pass three tests in the Peak Detection Module 
in order to qualify as a peakc^ tall- enough?, local'-maxima?, and shape—of -a-peak?. 

The tall-enough? Criterion 

In order to pass the t all - enough? criterion, the height of the potential peak must be 
sufficiently greater than the average height of the noise in the power spectrum. 

Although windowing the input signal helps minimize spectral leakage and adjust¬ 
ing the range of the HP 5182a waveform recorder helps reduce the quantization error, 
there still exists a nontrivial amount of noise scattered throughout the power spectrum. 
Therefore, to prevent erroneous peak detection, it is important that the tall-enough? 
module verify that the average height of each potential peak is sufficiently greater than 
the average height of the low-level noise. To determine what is “sufficiently greater,” the 
peak-detector module has a variable called fuzz-threshold. Any point with height less 
than this fuzz-threshold will not qualify as a peak point. The value of fuzz-threshold 
is a function of quantization error, arithmetic error and spectral leakage. 

The local-maxima? Criterion 

In addition to being taller than the noise floor, the potential peak must also be a local 
maximum. This means that the points to the right and to the left must be shorter than 
the potential peak. If, for example, the point to the right is taller, then the program 

2 The approach I’ve used for peak detection is very similar, though not based upon, the method used 
by Sreenivas and Rao[SR79] for pitch extraction. 
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discards the original point and adopts this right-hand point as a new peak candidate. It 
then runs the same tests on this new point to see if it meets all the criteria for a peak. 

The shape-of-a-peak? Criterion 

Besides being a local maximum and sufficiently tall, the potential peak must also have 
the shape of a peak. An ideal peak is characterized by a sharply monotonically-increasing 
slope and a sharply monotonically-falling slope. 

The sharp-enough? module compares the ratio of the potential peak’s height to that 
of a neighboring point. This ratio essentially gives the slope of the peak. If the slope is 
greater than a certain threshold value, then the peak is sharp-enough?. The advantage 
of using this method is its simplicity. The disadvantage is that the ratio is fairly sensitive 
to variations in the neighboring value. For example, it is possible to get a large ratio if 
the potential peak value is large or if the neighboring value is small. The first is desirable 
whereas the second is not. Yet the program avoids erroneous peak detection by imposing 
the tall-enough? condition, since the height of the potential peak must be greater than 
the height of the low-level noise. This prevents the peak-detector module from picking 
up extraneous short peaks merely because the neighboring value is smaller than usual. 

Periodogram Averaging or Bartlett Estimator 

For the many test cases I have run, the power spectrum was dean enough so the detector 
module could detect at least the salient peaks in the spectrum. However, in the event 
that the power spectrum, which is normally computed by taking a single FFT of the en¬ 
tire input signal, is too noisy and peak detection becomes too difficult, I have written an 
additional module capable of producing a cleaner power spectrum. Althou gh computa¬ 
tionally more expensive, this optional module can produce a dearer power spectrum using 
the periodogram averaging or Bartlett’s estimator technique[OS75, AN83, SS75, RG75] 
as outlined below. 

To get cleaner looking power spectra, 


1) take signal | 
of 16384 pts. 
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-1st chunk—| 

I-2nd chunk-1 

I—3rd chunk—| 

etc. 

3) FFT the chunks 

4) Average the results. 

More formally, in the periodogram averaging approach, a data sequence x(n), 0 < 
n ^ N ~ is divided into K segments of M samples each so that N = KM; i.e., we 
form the segments 

x(‘>(n) = x(n + iM - M), 0 < n < M - 1, 1 < * < (3.2) 

and compute the K periodograms 

i M ~ l 2 

/(0m (^) = , 1 < i < K (3.3) 

nso 

The spectrum estimate is then defined as 

1 K 

Bxx(u) — y^fj/(k>) (3.4) 

t=l 

Choosing the optimal chunk size is important to ensuring a clean power spectrum. In 
general, it is not obvious how to choose the optimal chunk size. Increasing or decreasing 
the chunk size introduces tradeoffs which must be weighed. In general this choice may 
be guided by prior knowledge of the signal under consideration. Refer to [OS75, RG75, 
OS75, AN83]. For the purposes of this program, an overlap of 50% works best[OS75]. 

3.2.2 Finding the precise location of a peak 

Once the program has identified a peak, it must record the peak’s location so the classi¬ 
fication component can determine its relationship to other peaks. 


2) divide signal |— 
into overlapping 
chunks of size 
2000 points 
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To compute the Frequency of a peak, the program uses Equation 3.5, 


( Frequency \ _ /_ Index _\ 

\Frequency Span) \ Length of Power Spectrum) \ • ) 

where Frequency Span and Lengthof Power Spectrum are given, and Index is computed 

by the peak-detector module. If each data point is in a one-to-one correspondence with 

each frequency point in the power spectrum, then equation 3.5 provides a simple and 

efficient way to determine the location of a peak. However, in some cases, there is 

far less than a one-to-one correspondence between points and frequencies. The forced 

negative-resistance oscillator exhibits most interesting behavior at frequencies from 50 

khz to 500 khz. Since the power spectrum in this program can be represented by at most 

8192 data points, 3 and the power spectrum must span 500 khz in order to reveal the 

interesting behavior, there is a 1 : = 1 : 6.103515625 correspondence between 

each data point and frequency. In other words, each data point in the power spectrum 

represents 6.103515625 hertz. As a result, the program might, in the worst case, be off 

by 6.103515625 hertz when it returns the location of a peak using equation 3.5. 

To alleviate the problem, I have implemented an additional module that computes 
more precisely the location of peaks based on the initial, rough estimate. 4 The strategy 
is as follows. 

First, for each peak, magnify the power spectrum around the peak using the chirp 
^-transform algorithm. See figure 3.6. The advantage of the chirp 2 -transform as it 
is implemented in the 6.003 software is that it allows the user to select the starting 
and ending frequencies, the incremental step size and the desired number of points for 
the magnified spectrum, unlike the fast Fourier transform procedure which does not 
afford such flexibility. In this program, I used 2048 data points to represent a magnified 

3 The HP5182& Waveform Recorder/Generator can capture a maximum of 16384 points of a single 
signal. Although the 6.003 FFT procedure outputs a spectrum of 16384 data points, half of these points 
provide redundant information since the power spectrum for the negative frequencies is merely a mirror 
image of the spectrum for the positive frequencies. Thus, after discarding the negative frequencies, we 
are left with 16384/2 = 8192 data points in the power spectrum. 

4 Finding the precise location of a peak is extremely important for classification purposes because of 
the integral relationships between subharmonic frequencies. Section two of this chapter will discuss this 
issue in greater depth. 
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frequency span of 30 hertz. This means that each peak location can be off by at most 
30/2048 = 0.0146484375 hertz, a significant improvement over 6.103515625. 

Second, take the tallest three points of the magnified spectrum, and fit them to a 
parabolic curve. Then use parabolic interpolation to find the center frequency of the 
peak. 5 In general, interpolation is necessary because the peak’s center may lie between 
two data points, and the use of a parabola is effective because peaks, when magnified, 
have shapes similar to inverted parabolas. See figure 3.6. 

More specifically, the equation for a parabola is: 


Ax 2 + Bx + C = y 

(3.6) 

Fitting the three points, (£ 1 , 2 / 1 ), (® 2 ,Sfe)> and (£ 3 , 2 / 3 ), to a parabolic 
solving three simultaneous equations in three unknowns. 

curve requires 

Ax\ + Bxi + C = yi 

(3.7) 

Ax\ + Bx-i + C = 1/2 

(3.8) 

Ax § + Bx 3 + (7 = 2/3 

(3.9) 

Solving for A, B, and (7, we get 


p - (jfe - yi)(xl - x\) - ( 2/2 - 2 /i )(*3 - * 1 ) 

(£3 - £i)(x^ - £j) - (£2 - £l)(>3 “ £ 1 ) 

(3.10) 

, {y* - y\) - b(x 2 - x x ) 

(*2 - *?) 

(3.11) 

(7 = 2/1 — -<4£i - Bx 1 

(3.12) 


5 This technique was suggested by Gerry Sussman. 
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f [khz] 


Figure 3.6: Magnification of Peak, (a) The top plot illustrates the full spectrum on a 
linear scale, (b) The bottom plot displays a magnified version of (a) around the peak. 

Once the values for A , B, and C are known, then the center frequency can be read¬ 
ily computed. This center frequency is nothing more than the maximum value of the 
parabola obtained by taking the derivative with respect to x of equation 3 . 6 . 

2Ax + B = 0 (3.13) 

Solving for x, 

JJ 

Center Frequency =-- ( 3 . 14 ) 

The third and final step is to take this center frequency of the peak and rescale it 
so it has meaning relative to the entire frequency span, not just the magnified span. 
Figure 3.6 and figure 3.7 illustrate the process. 

Magnifying each peak and interpolating to find the center frequency proves very 
effective. The result is a peak location much more accurate than that obtained using 
equation 3.5. The tradeoff, however, is that the chirp z-transform is computationally 
expensive. If there are many peaks in the power spectrum, then determining the precise 
location of all the peaks becomes very slow. 
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Figure 3.7: Parabola Fitting and Finding Center Frequency. 

3.3 Classification Component 

In the previous section, I describe the traits the peak-detection component looks for 
when identifying peaks and how it computes their locations in the power spectrum. 
Once the peak detector has gathered all the peaks and the frequencies at which they 
occur, it passes this list of peaks onto the classifier. In this section, I discuss the classi¬ 
fication component. The purpose of this module is to determine if the response exhibits 
harmonic, subharmonic, quasiperiodic or chaotic behavior. I begin this section with a 
precise definition of the terms. Following the definitions is an explanation of each of the 
modules within the classification component. 

3.3.1 Defining the Terms 

When looking at the power spectrum, the program will use the following simple defini¬ 
tions to identify the system’s behavior. Refer to figure 3.8 and figure 3.9 for examples. 

• Harmonic response — peaks in the power spectrum may occur only at the drive 
frequency, /drive > and integer multiples of the drive frequency. 

• Subharmonic response of order n — peaks in the power spectrum may occur only 
at £ • /drive and integer multiples of £ • /drive- 

• Quasiperiodic response of degree n — peaks in the power spectrum may occur 
only at n fundamental frequencies and integer combinations of these fundamental 
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frequencies. For example, a quasiperiodic response of degree 2 may have peaks 
only at 

frequencies = a • f\ + b • / 2 

where f\ and / 2 are not rational multiples of'each other and a and b are integers. 6 

• Chaotic response — the power spectrum has a lot of broadband noise in addition 
to some sharp peaks. The level of noise in the chaotic case is much greater than 
the level of noise for the harmonic, subharmonic, and quasiperiodic cases. 


3.3.2 Chaos Module 

The Classification Component of the program is made up of several modules, the first of 
which is the chaos module. Refer to Figure 3.10. The purpose of the chaos module is to 
identify chaotic behavior using evidence from the power spectrum. 

As mentioned in the definition above, the power spectrum of a chaotic response is 
characterized by broadband noise in addition to some sharp peaks. In addition, the 
level of noise is much greater in the chaotic case than the low-level noise observed in 
the periodic and quasiperiodic responses. Figure 3.11 provides an example of a power 
spectrum for a chaotic response. 

To recognize chaotic behavior in a power spectrum, I tried computing the average 
height of the spectrum and comparing this with the average height of a known, non- 
chaotic power spectrum. If the average height of this potentially chaotic response is 
significantly greater than that of a known, nonchaotic response, then I classified the 
system as chaotic. However, this proved to be an unreliable test for chaos. 

An alternative strategy for discerning chaos in the power spectrum involves searching 
for sharp, distinct peaks in the spectrum. If the peak detector returns no sharply defined 
peaks, we may conclude that the spectrum contains only broadband noise and hence 
classify the behavior as CHAOTIC. However, the fact that a response has sharply defined 
peaks in its spectrum does not preclude it from being classified as chaotic. As indicated 
in the definition of chaos, it is perfectly possible for a chaotic response to have some 

®The task of identifying quasiperiodicity is simplified slightly by the fact that the degree of a quasiperi¬ 
odic response seldom exceeds 2. According to Thompson and Stewart[TS86], page 196, “A remark¬ 
able theorem in qualitative dynamics, proved recently by Ruelle, Takens, and Newhonse, suggests that 
quasiperiodicity of degree greater than 2 may be difficult to observe in nature.” 
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Figure 3.9: Examples of power spectra on a logarithmic scale: (a) a harmonic response 
with fdrive = 40000 hz; (b) a second-order subharmonic response with /drive = 60000 
hz; (c) a quasiperiodic response with /drive — 69000 hz; and (d) a chaotic response with 
/drive = 66000 hz. 
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Classification 



Figure 3.11: Power spectrum of a chaotic response. 













0 khz 25 khz 50 khz 75 khz 100 khz 1 

Figure 3.12: Power spectrum of a second-order subharmonic response with /drive = 50 
khz. 

sharp peaks. I have found that the power spectrum does not easily provide conclusive 
information regarding the chaotic versus non-chaotic nature of a signal. 7 In the next 
chapter, I present the autocorrelation technique as a promising alternative approach to 
distinguishing between quasiperiodicity and chaos. 


Referring back to figure 3.10 again, we find that if the program has concluded the 
response is not chaotic, its next step is to distinguish between quasiperiodic and periodic 
behavior. Before discussing how the program makes this distinction, I wish to first 
explain how the program determines the order subharmonic of a periodic response. An 
understanding of how the program computes subharmonics will aid in the understanding 
of how it distinguishes between periodicity and quasiperiodicity. 


3.3.3 Subharmonic Module 

Assuming we already know the response is periodic, the next goal is to determine its 
order. Figure 3.12 shows the power spectrum of a second-order subharmonic response 
whose drive frequency is 50 khz. 

We recognize this response as second-order subharmonic because peaks occur only at 
integer multiples of: 

As it turns out, the problem of distinguishing between quasiperiodic and chaotic behavior is fairly 
difficult[UA81, KKC]. 


[khz] 
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(3.15) 


- • /drive where n = 2 and f drive = 50 khz 

71 

at frequencies 

25 khz — ^ • fdrivci 
50 khz — ^ • fdrive > 

75 khz — ^ • fdrive > and 
100 khz = | • / rfr ,„ e 

Information contained in the first 50 khz of the power spectrum is sufficient for us to 
classify the system’s response as second-order subharmonic. However, the information 
contained in the second segment from 51 khz to 100 khz also provides valuable supporting 
evidence for our hypothesis. In particular, the peak at 75 khz confirms the information 
contained in the 25 khz peak, while the peak at 100 khz confirms the information con¬ 
tained in the 50 khz peak. Scientists, especially when confronted with a noisy power 
spectrum, often search for this kind of redundant evidence to confirm or disconfirm their 
hypothesis. To take advantage of this duplicate information, I have chosen a multi-set 
representation for the list of peaks. This multi-set representation lumps together si mil ar 
peaks much as a scientist looking at a power spectrum would. For the power spectr um 
shown in figure 3.12, the list of peaks is 

((25000 . 10.4206) (50000 . 5.954) (75000 . 1.9) (100000 . 3.721)) 

I I_height of peak 

I 

location of peak 

where the first element of each set indicates the location of the peak (in hertz) and the 
second element indicates the height of the peak. This list of peaks will get transformed 
into the multi-set 
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((25000 . 2) (50000 . 2)). 

I I 

normalized number of occurrences 
location of 
peak 

with the first element of each set specifying the normalized 8 location of analogous peaks 
and the second element specifying the number of occurrences of such analogous peaks. 
Thus, in this case, the (25000 . 2) entry indicates that there were two occurrences 

of the 25000 hertz “type” peak in the power spectrum, one at 25000 hertz and one at 
75000 hertz. Similarly, there were two occurrences of the 50000 hertz “type” peak, one 
at 50000 hertz and one at 100000 hertz. 

Once the program has computed this multi-set of remainders, it then transforms this 
multi-set of remainders into a multi-set of hypotheses by dividing the drive frequency 
by each normalized location. Upon performing this division, we get a new multi-set of 
hypotheses. 

((2 . 2 ) (1 . 2 )) 

I I_number of occurrences 

I 

hypothesized 
order subharmonic 

To compute the final order of the system’s response, the program takes the least com¬ 
mon multiple of all the hypothesized orders. In this example, (1cm 2 1) » 2, indicating 
a second order subharmonic response, which is correct. 

If the peak detector fails to return the precise location of the peaks, but instead 
returns the following list: 


*The normalized location is computed by taking the (peak location ) sod (drive frequency). 


((25000.003 .2.8) (50000.012 . 10.206) (75000.001 .3.1) (100000 . 6.229)) 

I I_height of peak 

I 

location of peak 

the program can still take advantage of the multi-set representation. This is true because 
the program tries to determine what the true location of the peak should have been. More 
specifically, if the computed location of the peak were 25000.003, but the true location is 
25000, then the program tries to justify renaming the 25000.003 peak to a 25000 peak. 

If the computed peak at 25000.003 is dose-enough? to the ideal peak location at 
25000, then the program goes ahead and renames the 25000.003 peak to a 25000 peak. 
Computing the close-enough? threshold level is easy. Recall from section 3.2.2, we 
know how to determine the maximum deviation of a computed peak location from its 
true location based upon the frequency span and the number of points in the power 
spectrum. Thus, if a computed location is off from the ideal location by less than this 
maximum deviation, then the difference in distance between the ideal location and the 
computed location could logically be explained by the coarse granularity of the power 
spectrum. 

In the event that the peak detector erroneously picks up a stray spike, the correct clas¬ 
sification of the system’s response becomes very difficult. The multi-set representation, 
however, helps the program to distinguish between true peak and stray spike. 9 

For instance, if the peak-detector module accidentally detected an unusually tall stray 
spike at 847.457627119 hertz, then we might get the following for the list of peaks: 

((847.457627119 . 1.04) (25000 . 9) (50000 . 4) (75000 . 2.3) (100000 . 3)), 


’Windowing the input to reduce spectral leakage, selecting the appropriate input range on the 
HP5182a waveform recorder and utilizing Bartlett’s periodogram averaging technique discussed above 
reduce the low-level noise in the power spectrum which minimizes the probability of erroneous peak 
detection by the Peak Detector Module. But, the multi-set representation helps the program recognize 
stray peaks after the Peak Detector has already identified these strays as peaks. 
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with a multi-set of remainders of: 


((847.457627119 . 1) (25000 . 2) (50000 . 2)), 


and a multi-set of hypotheses of: 

((59 . 1) (2 . 2) (1 . 2)) 

Taking the least common multiple of the 59, 2 and 1, we get 118. At this point the pro¬ 
gram is smart enough to recognize that a 118th order subharmonic response is 1) nearly 
impossible to measure because of the imprecisions of the instruments, and 2) is highly 
unlikely even if the instruments were extremely precise. Instead of returning a conclusion 
of (SUBHARMONIC 118), it goes back to the multi-set of remainders and notices that there 
was only a single occurrence of (847.457627119 . 1.04). Since there was only a single 

occurrence of the 847.457627119 “type” peak, this suggests that the 847.457627119 peak 
may be a stray. (Perhaps the peak-detector module erroneously picked it up.) The pro¬ 
gram then knows to toss out this single, stray peak and recomputes the order until the 
answer seems reasonable. In this case, throwing out the (847.457627119 . 1) term 

from the multi-set of remainders allows the program to recognize the system’s response 
as second-order subharmonic, as it should. 

Quasiperiodic-or-Periodic Module 

Now that we know how the subharmonic module classifies the various order subharmonic 
responses, we can begin to discuss how it distinguishes between quasiperiodicity and 
periodicity. 

Recall from section 3.3.3 that if there was a stray peak in the power spectrum, which 
caused the computed order to be extremely high, the program was smart enough to throw 
out this peak and recompute the order until it arrived at a reasonable order. For periodic 
responses there should be no more than a few of these stray peaks. In the section on 
the subharmonic module we saw how nicely the 25 khz peak got lumped with the 75 khz 
peak and how the 50 khz peak got lumped with the 100 khz peak when we normalized 
with the drive frequency of 50 khz, leaving no stray peaks. 
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Figure 3.13: Power spectrum of a quasiperiodic response with two fundamental frequen¬ 
cies at u >3 and u>i. 

The power spectrum of a quasiperiodic waveform, with sharp peaks at multiple fun¬ 
damental frequencies and integer combinations of these fundamental frequencies, will be 
much more likely to have a greater number of single occurrence peaks than a periodic re¬ 
sponse. Referring to the quasiperiodic power spectrum in figure 3.13[TS86], we see peaks 
at the two fundamental frequencies, u>i and u> 3 , along with peaks at some of their integer 
combinations, 2u$ — — u 3 , and — o; 3 . Normalizing all the peaks at %j$ — u>i, 

u>i — u> 3 , W 3 , %}\ — u >3 by the drive frequency yields a different answer almost every time 
because each peak is not an integer multiple of another peak. Thus, it becomes much 
more difficult to lump together analogous peaks as we could in the periodic case when 
peaks tended to be multiples of another rather than sums and differences of one another. 

To distinguish between quasiperiodicity and periodicity, the quasiperiodic-or-periodic 
module notices how many single occurrence peaks occur in the multi-set of remainders. If 
more than a few single occurrence peaks exist in a spectrum , 10 this suggests the system’s 
response may be quasiperiodic . 11 

10 For the purposes of this implementation of the program, “a few” means two single occurrence peaks. 
This is a fairly tight constraint. However, based upon the examples seen so far, the power spectra of 
periodic signals tend to have sharp peaks only at integer multiples of each other, producing very few 
single occurrence peaks in the multi-set. 

11 Actually, it would be nice if we could easily group together analogous peaks in the quasiperiodic 
spectrum just as we can in aperiodic spectrum. For example, the peaks at 2 u >3 -wi.wi - 013 , 0 *, wi, and 


If the quasiperiodic-or-periodic module deduces that the power spectrum indi¬ 
cates quasiperiodic behavior, the program returns (QUASIPERIODIC); otherwise, it pro¬ 
ceeds to the subharmonic module and computes the order subharmonic as described in 
the previous section. 


2u>i — M3 should all be lumped together since they are all integer combinations of the drive frequency, 
Mi, and the other fundamental frequency, M3. The problem is that in order to recognize the analogous 
peaks, we must first be able to recognize all the fundamental peaks, which proves difficult. Even with 
the two fundamental frequencies known, there may be too many combinations to detect. 
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Chapter 4 


The Time-signal Analyzer 

In the previous chapter, we discussed how the power-spectrum data can be used to 
classify a system’s response. In most cases spectral analysis can accurately categorize 
the behavior of a waveform. However, as mentioned in chapter 2, the power spectrum 
provides only one perspective of the system’s behavior. When the noise level in the spec¬ 
trum increases, peak detection and behavior classification become more difficult, thereby 
increasing the likelihood of erroneous classification. To supplement the information pro¬ 
vided by the power-spectrum data, the program also looks at the time-signal data. The 
purpose of this chapter is to present the time-signal analyzer, which identifies the circuit’s 
behavior using time signal information. 


4.1 Time-signal definitions 

Before discussing the details of implementation, I first present the definitions of harmonic, 
subharmonic, quasiperiodic, and chaotic from a time-domain perspective. 1 These are the 
definitions used by the time-signal analyzer to categorize the behavior of a signal. 

• Harmonic response — a periodic response with the same frequency as the driving 
frequency. T drive — Tresponaet where T — period,. 

• Subharmonic response of order n — a periodic response whose frequency is some 
fraction of the drive frequency. Or T reaponae = n • T drive , where T = period. 


The previous definitions given in chapter 3 were from a power spectrum perspective. 
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• Quasiperiodic response a response that is the sum of a couple of periodic signals, 
whose periods are not rational multiples of each other. The signal is not periodic, 
but rather, almost periodic. 

• Chaotic response — a response that cannot be simply be characterized. The time 
signal does not repeat itself in a finite amount of time. 

4.2 Basic strategy of the time-signal analyzer 

This section describes the strategy used by the time-signal analyzer to identify the be¬ 
havior of a recorded waveform. 

4.2.1 Strobe technique to determine order subharmonic 

The time-signal analyzer first assumes the recorded waveform is periodic. It then deter¬ 
mines the order subharmonic by strobing the signal at integer multiples of the drive period 
and recording the value of each strobe. Since the program knows the drive frequency, 2 
it can compute the drive period using the relation T drive = l/f drive . If all strobed values 
are within a given epsilon of each other, then the signal is said to be harmonic, since 
it repeats itself at intervals of a period. If every other strobe point has the same value, 
then the waveform is second-order subharmonic. If every third strobe point has the same 
value, then the waveform is third-order subharmonic, etc. Figure 4.1 through figure 4.4 
illustrate the notion of strobing the signal to determine order subharmonic. In each of 
the figures, the top graph displays the forcing sinusoid, and the bottom graph shows the 
response. 

In the actual implementation, I examine ten values spaced at multiples of the drive 
period. If each of the later nine points are within a small distance, epsilon, of the first 
point, then the signal is periodic, more specifically harmonic. If not, the program looks 
at ten values 3 spaced at twice the drive period. Once again if the later nine points are all 

Recall the inputs to the entire analyzer program are the ranges of frequency and voltage over which 
to run the analysis. 

For greater assurance of accuracy, you can increase the number of strobed points. Increasing the 
number of strobed values from ten to say fifty means that fifty points must be “dose enough” to one 
another as opposed to ten. This is a much tighter constraint. However, if all fifty points are “close 
enough,” then you can be more certain that the deduced subharmonic order is correct. 
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Figure 4.1: Strobing a harmonic response. Notice how the response waveform repeats 
itself at every strobe. 



t [sec] 


Figure 4.2: Strobing a second-order subharmonic response. Notice how the response 
waveform repeats itself once for every two times the driving sinusoid repeats itself. 
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'*■ t [sec] 





Figure 4.3: Strobing a quasiperiodic response. 


t [sec] 


t [sec] 



Figure 4.4: Strobing a chaotic response. Notice how the response never repeats itself, 
even with an infinite number of strobes. 
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within a short distance of the first point, then the program calls this signal subharmonic 
of order two, and so on up to some large order subharmonic, like thirty-three, 4 If the 
order subharmonic is greater than thirty-three, then the behavior of the response is 
probably quasiperiodic or chaotic, and the program classifies it as such. 

To be truly precise, the time-signal analyzer should probably compare each strobed 
value with every other strobed value, no, jus, the first one. In my implementation, 
the determination of the order subharmonic is highly dependent upon the value of the 
first strobed point. If for some reason the value of this first point has a large variance 
m comparison to the other nine points, then the analysis returned by the time-signal 
analyzer may be inaccurate. However, the advantage of this implementation is that i, 
.« simpler and more efficient. Moreover, it is unclear how much certainty we gain by 
checking extra point values. For the examples I have run, there doe, no, seem to be 
a sacrifice in the accuracy of the conclusions resulting from comparisons with the later 
nine points with the firs, point only. If, however, this does prove to be a problem, we 
can easily modify the time-signal analyzer so it performs a point-by-poiut checking. 

Finally we must address the issue of how the program computes epsilon. Recall 
that epsUon is the threshold value below which two strobed points are considered close 
enough and above which two strobed points are not eke, enough. The time-signal 
analyzer computes the epsilon value by taking a certain fraction of the signal span. More 
specifically, the program computes the difference between the maximum and minimum 
values in the rime signal. This is the spun of the signal. I, then multiplies the span by 
a certain fraction, in this case 0.10. Wha, this means is if a subsequent strobed point 
vanes from the first strobed point by no more than 10% of the spun, then the two points 

are essentially close enough, and the signal has a chance of repeating itself in a periodic 
or quasiperiodic manner. 

high.,,h» tUrtyim .t Tf fol *“■ “ *■“* subhstmoiiics 

imtmmsuw and “ ‘ h ' 
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Figure 4.5: Autocorrelation of a Quasiperiodic Response. 

4.2.2 Autocorrelation technique to distinguish between quasiperiodic¬ 
ity and chaos 

The previous section described how the time-signal analyzer computes the order sub¬ 
harmonic of a periodic response. But how does the analyzer recognize a chaotic or 
quasiperiodic response? One possible technique for distinguishing between a quasiperi¬ 
odic and a chaotic response is to compute the autocorrelation of the recorded waveform. 5 
If the signal remains fairly highly correlated with itself through time, it is periodic or 
quasiperiodic. See Figure 4.5[Lin88]. If, on the other hand, the correlation values are 
very low as shown in figure 4.6[Lin88], then the signal is chaotic. 

An autocorrelation value of zero is said to be uncorrelated , while a large autocorrelar 
tion value indicates high correlation. When we look at figure 4.6, it makes sense that the 
chaotic signal, which never repeats itself, has low correlation values for large time delays 
and high correlation when the time delay equals zero, since the signal correlates exactly 
with itself, but not very well with the shifted versions of itself. The quasiperiodic signal 
has the property that it is “almost periodic”; that is, the signal almost repeats itself. 
Therefore, as shown in figure 4.5, its autocorrelation function has fairly large values in 

This strategy of using the autocorrelation to distinguish between quasiperiodicity and chaos was 
suggested during a discussion with Professor Paul Linsay of the MIT Physics Department. 


59 





1.0 



0 4 8 12 16 

Time Delay (10 3 drive periods) 


Figure 4.6: Autocorrelation of a Chaotic Response. 

comparison to the the chaotic case, but as the time delay grows larger, the amount of 
correlation also decreases. 

The notion of using the autocorrelation of a signal to distinguish between quasiperi- 
odic and chaotic responses looks promising. However, I have not yet implemented this 
component yet. Paul Linsay[Lin88], in his studies of quasiperiodicity and chaos, has used 
the autocorrelation technique extensively to distinguish between the two types of behav¬ 
ior and finds this method highly effective. The implementation of the autocorrelation 
component of the time-signal analyzer has been left for future work. 

4.3 The Resolver 

What happens if the conclusion produced by the power-spectrum analyzer disagrees with 
the conclusion of the time-signal analyzer? In such cases, the resolver module must re¬ 
solve the conflict. To decide which one of the two conclusions to accept, the resolver 
will look at contextual information. If, for example, the power-spectrum analyzer iden¬ 
tified the behavior as CHAOTIC, and the time-signal analyzer classified the behavior as 
SUBHARMONIC ORDER 2, the resolver will look at the neighboring points to the east and 
to the south. If both are SUBHARMONIC ORDER 2, then the resolver will assume that the 
power spectrum was extraordinarily noisy, and it will reject the CHAOTIC analysis and 


60 


accept the SUBHARMONIC ORDER 2 analysis. This strategy is based upon the assump¬ 
tion that within a given local region, the behavior of a system will tend to behave in a 
similar manner. The validity of this assumption depends upon the step size or granu¬ 
larity with which we have chosen to run the program and the rapidity with which the 
circuit changes behavior. Contextual information provides just one method of resolving 
conflicting analyses. The final chapter of this thesis suggests others. 
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Chapter 5 


The High-level interpreter 


In the previous two chapters of this thesis, I discussed the classification module. Recall 
that the primary purpose of the classification module is to characterize the behavior of 
the nonlinear system for every point in the two-dimensional parameter space. In this 
chapter, I present the high-level interpreter. Its purpose is to take the detailed parameter- 
space graph produced by the classification module and to provide a high-level, textual 
explanation describing any interesting, qualitative properties of the circuit’s behavior. 

For example in chapter 1, we saw that the high-level interpreter generated the following 
textual explanation given the parameter-space graph shown in figure 5.1. 

HIGH-LEVEL INTERPRETATION OF PARAMETER-SPACE DIAGRAM 

The parameter space diagram characterizes the behavior of the 
forced negative-resistance oscillator as the frequency and amplitude 
of the driving sinusoid are varied from 40000 to 70000 hz and from 
1.5 to 6.5 volts respectively. For amplitudes between 1.5 and 6.5 
volts and frequencies between 40000 and 56000 hertz, the system ex¬ 
hibits predominantly first order subharmonic behavior, while for 
amplitudes between 1.5 and 5.75 volts and frequencies between 41000 
and 68000 hertz, the system exhibits predominantly chaotic behavior, 
while for amplitudes between 1.5 and 6.5 and frequencies between 
55000 and 70000 hertz the system exhibits predominantly second order 
subharmonic behavior. Aside from these large areas of uniform behavior, 
the system’s response passes briefly through subharmonics of order 
3, 4, 5, 7 and 9 at isolated points throughout the parameter-space 
graph. 
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Figure 5.1: Sample parameter-space diagram. 
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The large region of first order subharmonic behavior is bounded 
by the bifurcation curves 1-2 and 1-C, where C is a region consisting 
of primarily chaotic behavior, but littered with sporadic occurrences 
of other kinds of behavior. Similarly, the large region of second 
order subharmonic behavior is bounded by the bifurcation curves 2-1 
and 2-C, while the region of chaos is bounded by C-l and C-2 bifur¬ 
cation curves. 

Finally, based upon the data gathered so far, the forced neg¬ 
ative-resistance oscillator shows no clear signs of period-doubling 
bifurcations as described by Feigenbaum and others. 

In order to provide such a high-level description of the circuit’s behavior, the inter¬ 
preter must identify large regions of uniform behavior as well as the small regions of 
sporadic behavior. It must notice the types of boundaries which border each region, and 
finally it must recognize such well-studied patterns as period-doubling routes to chaos 
and period-adding routes to chaos. How the program accomplishes all these tasks is the 
topic of this chapter. 

The first section in this chapter discusses how the high-level interpreter scans the pa¬ 
rameter space to produce the first paragraph of the textual explanation, which describes 
the different regions of behavior. The second section of this chapter explains how the 
interpreter uses computer vision techniques to detect the edges of a region and to locate 
the bifurcation curves that border the regions of uniform behavior. The third describes 
how the interpreter uses knowledge about general nonlinear dynamics to identify such 
typical patterns such as period-doubling cascades and period-adding cascades in the pa¬ 
rameter space. The fourth section talks discusses the symbolic description returned by 
the analyzer program and used as input to other programs. 


5.1 Growing regions of uniform behavior 

This section describes how the high-level interpreter identifies the large and small areas 
of uniform behavior in the parameter space. By utilizing the computer vision technique 
of growing regions of uniform behavior[Hor86], the interpreter can locate the regions of 
qualitatively distinct behavior. The strategy for growing such regions is as follows. 
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The program begins with any point in the parameter space, say the origin. This 
becomes the seed point. It then tries to grow the seed in all possible directions: north, 
south, east and west. If, for instance, the behavior of the point to the north is the same 
as the behavior of the seed point, then the area of this region expands to include the 
northern point. In other words, any point contiguous to the seed point and having the 
same behavior as the seed point becomes incorporated into the growing region. After 
this northern point has been CLAIMED by a given region, it can no longer be CLAIMED by 
another region. Instead, it will serve as a new seed point for further growth of the region. 
When the region cannot grow any further in any direction, the program systematically 
scans the parameter space point by point for the first UNCLAIMED point to serve as a new 
seed point for a new region of behavior. Finally, when all the points in the parameter 
space are CLAIMED, we have a parameter space in which all contiguous points having the 
same behavior are grouped together as desired. See Figure 5.2. Given this grouping of 
points, the high-level interpreter can easily generate statements similar to those shown 
in the first paragraph of the sample textual explanation. It can describe the large regions 
of uniform behavior, draw attention to isolated regions of sporadic behavior, make note 
of general changes in behavior as an input parameter is varied, and comment upon the 
frequencies and voltages at which these regions of behavior occur. 

5.2 Recognizing boundary types of a region 

Generating the textual explanation for the second paragraph, however, is not as straight¬ 
forward. Rather than requiring only local information about the behavior of neighboring 
points, the program now needs more global information. For instance, it must be able 
to look at the parameter-space diagram in figure 5.1 and recognize that the large region 
of first-order subharmonics is bounded by two kinds of borders: a 1-C border, and a 1-2 
border, while the region of second-order subharmonics is bounded by a 2-1 bifurcation 
curve and a 2-C bifurcation curve. The high-level interpreter must also recognize that 
the streaks of 3rd, 4th, 5th, and 9th-order subharmonics in the chaotic region do not 
change the predominance of chaotic behavior in the middle frequencies. In other words, 
the behavior in the middle region should be viewed as mostly chaotic with sporadic 
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Figure 5.2: Lumped regions of uniform behavior. 
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occurrences of other types of behavior, rather than as ten separate regions of behavior. 

Figure 5.3 illustrates the overall strategy of the High-Level Interpreter in recognizing 
the boundary types of a region. In general the task may be decomposed into two main 
parts. The first, which is illustrated in the first column of blocks in Figure 5.3, is to 
identify the large regions of uniform behavior. The second, which is illustrated in the 
second column of blocks in Figure 5.3, is to detect the bifurcation curves bordering each 
region of uniform behavior. 

5.2.1 Identifying large regions of uniform behavior 

In cases where the parameter-space diagram is not littered with small, isolated occur¬ 
rences of sporadic behavior, recognizing the regions of uniform behavior is easy. The 
program need only grow points of similar behavior using the algorithm described in the 
previous section. If, however, the parameter space is Uttered with brief occurrences of 
sporadic behavior or “blemishes,” the interpreter attempts to remove these local blem¬ 
ishes, since they distract from the overall picture. A blemish is defined as a small region 
in the parameter space having an area less than a certain fraction of the total param¬ 
eter space area. 1 To remove a blemish, the interpreter checks to see if the neighbor to 
the north belongs to a large region of uniform behavior. If it does, the interpreter will 
copy the behavior of the northern point into the behavior of the current point, thereby 
erasing the blemish. If the northern point does not belong to a large region of uniform 
behavior, the interpreter will check to the south, then to the east and then to the west. 2 
When all the blemishes in the entire parameter space have been removed, the program 
re-grows the regions of uniform behavior. After the second growth, there should be fewer 
distinct regions in the parameter space, but each region should have a greater area or a 
greater number of points in it. Figure 5.4 iUustrates the result of removing the blemishes 
from the parameter space in Figure 5.1. Notice how the program can now easily recog¬ 
nize the three large regions of uniform behavior: first-order subharmonic, second-order 

'Currently the threshold area is set at 3% of the total area. 

2 The problem of blemishes in the parameter-space graphs is similar to Horn’s problem of salt-and- 
pepper noise m a picture. To remedy the problem of “noise,” we both take advantage of contextual 
information. 
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Figure 5.4: Cleaned up parameter space without blemishes, 
subharmonic, and chaos. 

Identifying large regions of very messy behavior 

Now consider another parameter space diagram shown in figure 5.5 with a large region 
of first-order subharmonics centered about 5.25 volts and 46000 hertz, a large region of 
second-order subharmonics centered about 5.5 volts and 65000 hertz, a large region of 
chaos centered about 3 volts and 60000 hertz, and a large region of very messy behavior 
centered about 2.25 volts and 55000 hertz. 

The high-level interpreter should be able to categorize this last, messy region as a 
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Figure 5.5: Parameter-space diagram with region of very messy behavior near 2.25 volts 
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single region of very messy behavior, not as twenty-two tiny, distinct regions. Viewing 
these small regions as distinct would require the interpreter to compute the boundary 
types for each tiny region within this very messy area of the parameter space. Although 
this is possible, we would not gain very much in doing so. What the high-level interpreter 
does, then, is it converts all blemishes which are adjacent to another blemish into a 
new kind of region, “M”, which stands for “messy” region. Thus, the parameter-space 
diagram shown in figure 5.5 would get transformed into the parameter space shown in 
Figure 5.6. After all the blemish points have been converted to “M,” the Interpreter 
grows regions of uniform behavior once again. The result is the lumping together of 
points exhibiting qualitatively similar behavior. Notice how the messy region is now 
treated as a large, single region rather than many small regions. 

The justification for categorizing these blemishes into a single region is as follows. An 
investigator examining the behavior of this nonlinear, electrical circuit would notice that 
there is a large region of first-order subharmonics at the lower frequencies, a large region 
of second-order subharmonics at the higher frequencies, a large region of chaos in the 
middle frequencies, and another region of very messy behavior at the lower frequencies 
and amplitudes. The details of this last region are not as crucial as the fact that the 
general behavior in the region is messy. In addition, from the practical point of view, 
identifying the different boundary types for each tiny blemish can be a tedious process. 
Instead, it is preferable to identify, in some rough, approximate sense, the large regions 
of uniform behavior and then to compute the boundary type for each. 

5.2.2 Detecting bifurcation curves that border the regions of uniform 
behavior 

At this point, the high-level interpreter has finished identifying the large regions of 
uniform behavior. The next task is to detect, for each region of uniform behavior, the 
bifurcation curves that border the region. In other words, we want to be able to say 
the region of chaos in figure 5.1 is bounded on one side by a C-l bifurcation curve and 
on the other side by a C-2 bifurcation curve. To be able to make such a statement, the 
interpreter must first locate the edges of the regions in the parameter space. 

Once again I have leveraged off of a computer vision technique called binarization to 


71 



A [V] 







6.5V 

1 

1111111111111112222222222222222 




6.25V 

1 

1111111111111112222222222222222 




6.V 

1 

1111111111111111122222222222222 




5.75V 

1 

1111111111111111 

.22222222222222 




5.5V 

1 

1111111111111111 

. .2222222222222 

KEY: 



5.25V 

1 

111111111111111. 

. .2222222222222 




5.V 

1 

1111111111111... 

. .2222222222222 

1 = 

1st order 

subharmonic 

4.75V 

1 

1111111111111... 

.2222222222 

2 - 

2nd order 

subharmonic 

4.5V 

1 

111111111111.... 

.222222222 

M - 

messy behavior 

4.2SV 

1 

11111111111. 

.22222222 

. * 

chaotic 


4.V 

1 

1111111111M. 

.22222222 




3.7SV 

1 

lllllliillM. 

....22222222222 




3.5V 

1 

111111111MM. 

.222222 




3.25V 

1 

llllllilMMM. 

.222222 




3.V 

1 

1111111MMMM. 

.22222 




2.75V 

1 

111111MMMMM. 

.22222 




2.5V 

1 

11111MMMMMM. 

.2222 




2.25V 

1 

1111MMMMMMM. 





2.V 

1 

111MMMMMMMHM.... 

.222 




1.75V 

1 

11MMMMMMMMMM.... 





1.5V 

1 

1MMMMMMMMMMM.... 






- 

-1-1- 

-1-1- 

-> f [khz] 



40 50 

60 70 





Figure 5.6: Parameter Space Diagram with Messy Region Indicated by “M.” 
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Figure 5.7: Binarized parameter-space diagram for the region of first-order subharmonics. 

facilitate in edge detection[Hor86]. The strategy is as follows. Take a single region of 
uniform behavior and convert all points within this region to l’s and all other points to 
0’s. Next, scan the parameter space for 1 —* 0 or 0 —*■ 1 transitions. Any time you come 
across such a transition, record that point as an edge point. Also note the behavior of 
the adjoining point that caused this point to be an edge point. For example, referring 
again to figure 5.1, suppose we decide to binarize the region of first order subharmonics 
first. The result of binarizing the parameter space graph would produce a new graph as 
shown in figure 5.7. 

We next scan the space in figure 5.7 searching for 1 —*• 0 transitions. Scanning the 
top row, we encounter our first 1 —*■ 0 transition as the frequency is increased from 54000 
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hertz to 55000 hertz. Upon detecting this edge point, the interpreter records the point 
at 54000 hertz and 6.5 volts as an edge point and notices that the eastern point caused 
this point to be classified as an edge point. The program then looks up the behavior 
of this eastern point by referring to the original, unbinarized parameter space graph, 
notices it is a second-order subharmonic and records the behavior along with the edge 
point. Noting the behavior of the boundary-causing point is important because it allows 
the high-level interpreter to identify the boundary type. 

In general, binarizing each region of uniform behavior has the advantage of simplifying 
the process of edge detection. The high-level interpreter need only search for 1 —► 0 or 
0 —* 1 transitions. Once all the boundary points and the behaviors of the edge-causing 
neighbors have been recorded, we are left with a data structure containing a summary 
of all the boundary types of a region. Figure 5.8 illustrates the contents of a sample 
boundary data structure. 

The particular example given in figure 5.8 came from the parameter-space graph 
shown in figure 5.1. Notice how the data structure accurately describes the behavior in 
the the parameter space. It says the region of first-order subharmonic, labelled “Rl,” is 
bounded by three types of regions: chaotic, messy, and second-order subharmonic. The 
region of messy behavior, labelled “R2,” is bounded by two types of regions: chaotic, and 
first-order subharmonic. The region of chaotic behavior, labelled “R3,” is bounded by 
three types of regions: messy, second-order subharmonic, and first-order subharmonic. 
And, the region of second-order subharmonic, labelled “R4,” is bounded by two types 
of regions: chaotic and first-order subharmonic. Once the high-level interpreter has 
this data structure, generating the second paragraph of the textual explanation, which 
describes the bifurcation curves that border each region of uniform behavior, becomes 
straightforward. 

5.3 Recognizing Typical Patterns in the Parameter Space 

Finally, the purpose of the last paragraph of the textual explanation is to recognize 
and describe any interesting patterns observed in the parameter space, such as period¬ 
doubling cascades or period-adding cascades. As shown in figure 5.3, the two inputs 
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Specific Example: 


(((SUBHARMONIC 1 "Rl") 

((MESSY "R2") (CHAOTIC "R3") (SUBHARMONIC 2 "R4"))) 
((MESSY "R2") 

((CHAOTIC "R3") (SUBHARMONIC 1 "Rl"))) 

((CHAOTIC "R3") 

((MESSY "R2") (SUBHARMONIC 2 "R4") (SUBHARMONIC 1 "Rl"))) 
((SUBHARMONIC 2 "R4") 

((CHAOTIC "R3") (SUBHARMONIC 1 "Rl")))) 


General Format: 


( <region-A> 

( <bordering region Al> <bordering region A2>... <bordering region AN>) 
<region-B> 

( <bordering region Bl> <bordering region B2>... <bordering region BN>) 


<region-Z> 

( <bordering region Zl> <bordering region Z2>... <bordering region ZN>)) 


Figure 5.8: Boundary Representation Data Structure returned by High-Level Interpreter. 
Notice how the regions are numbered “Rl” to “RN” to distinguish among multiple regions 
with the same kind of behavior. 



to this module are: 1) the boundary representation data structure computed using the 
method outlined in the previous section, and 2) general knowledge about nonlinear dy¬ 
namics. Given these two pieces of information, the high-level interpreter has enough 
information to identify and explain any interesting patterns. As an example, the pro¬ 
gram searches for a period doubling by tracing through the boundary data structure 
searching for the pattern 1—►2—>4—* 8 —► 16 *■ chaos. In reality, the inter¬ 

preter does not actually trace through the full 1 -► 2 -+ 4 -+ 8 -► 16 -*• ... -► chaos 
sequence. Instead, if it can successfully trace through four contiguous regions where the 
subharmonic order of the current region is twice the order of the previous region, then 
the interpreter declares the existence of a period-doubling cascade. 

To expand the knowledge base of the interpreter, simply add new procedures to 
the high-level interpreter module. Each procedure must, of course, search out a given 
pattern such asl—>2—>-4—>8—► 16—chaos for the period-doubling route 
to chaos, l->2->3-^4->5-> chaos for the period-adding route to chaos, 

and periodic —► chaotic —► periodic —► chaotic —*• periodic —► chaotic —*• periodic ... for 

the alternating periodic-chaotic sequence. A particularly interesting route to chaos is 
called intermittency. To recognize an intermittent sequence to chaos, see if the system 
undergoes a transition from periodic behavior to chaotic behavior with occasional bursts 
of noise. Initially there should be long intervals of periodic behavior between short bursts 
of noise, but with increasing time the intervals between the bursts decrease; it becomes 
more and more difficult and finally impossible to recognize the regular oscillations of the 
periodic state[TS86]. 

5.4 Symbolic description 

Finally, after producing the parameter-space diagram and the accompanying high-level 
text explanation, the program returns a symbolic summary of the circuit’s behavior. Re¬ 
call from chapter 1 that the symbolic descriptor may be passed on to another procedure 
for further analysis of the circuit’s behavior. The other program might be interested in 
searching for more complicated patterns such as saddle-node bifurcations, flip bifurca¬ 
tions, or Hopf bifurcations[TS86]. Or the other program may be interested in locating 
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Specific Example: 


( t(((CHAOTIC) ;; first region of uniform behavior 

(((VOLTAGE 1) (DRIVE-FREQUENCY 57000)) ;; pts. in chaotic region 
((VOLTAGE 4) (DRIVE-FREQUENCY 56000)) 


((VOLTAGE 1) (DRIVE-FREQUENCY 50000)))) 

((SUBHARMONIC 2) ;; second region of uniform behavior 

(((VOLTAGE 2) (DRIVE-FREQUENCY 68000)) ;; pts. in subharm 2 region 

((VOLTAGE 2) (DRIVE-FREQUENCY 69000)) 

((VOLTAGE 1) (DRIVE-FREqUENCY 69000)))) 

((SUBHARMONIC 5) (((VOLTAGE 3) (DRIVE-FREQUENCY 56000)))) 
((SUBHARMONIC 7) (((VOLTAGE 3) (DRIVE-FREQUENCY 59000)))) 


((SUBHARMONIC 1) 

(((VOLTAGE 5) (DRIVE-FREqUENCY 51000)) 

((VOLTAGE 5) (DRIVE-FREQUENCY 50000)) 

((VOLTAGE 4) (DRIVE-FREQUENCY 50000)))) 

((SUBHARMONIC 9) (((VOLTAGE 4) (DRIVE-FREqUENCY 59000)))) ) 

(((SUBHARMONIC 1 "Rl") ;; boundary representation 

((MESSY "R2") (CHAOTIC "R3") (SUBHARMONIC 2 "R4"))) 

((MESSY "R2") 

((CHAOTIC "R3") (SUBHARMONIC 1 "Rl"))) 

((CHAOTIC "R3") 

((MESSY "R2") (SUBHARMONIC 2 "R4") (SUBHARMONIC 1 "Rl"))) 
((SUBHARMONIC 2 "R4") 

((CHAOTIC "R3") (SUBHARMONIC 1 "Rl")))) ) 


Figure 5.9: Sample symbolic description returned by the program. The first part of the 
symbolic descriptor describes the lumped regions of uniform behavior, while the second 
part provides information concerning the boundary representation. 
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Chapter 6 


Conclusion 


This thesis describes the design and implementation of a system that automatically 
measures and classifies the behavior of driven, nonlinear circuits and generates high- 
level, qualitative interpretations of their behaviors. In the concluding chapter, I would 
like to focus on the contributions and future work of this thesis. 

6.1 Contributions of Thesis 

6.1.1 Valuable Tool for Experimental Dynamicists and Other Investi¬ 
gators 

At present the analyzer program is still at an early stage of development. It has only been 
tested on a few sample circuits — the forced negative-resistance oscillator circuit[UA81], 
the series RLC circuit, and the devil’s staircase drcuit[KKC]. However, the results look 
promising. The analyzer program, along with the signal-processing module, the commu¬ 
nication module, the measurement instruments, and the observation equipment provide 
a valuable tool for experimental dynamicists and other scientists studying the behav¬ 
ior of complex, nonlinear, electrical circuits. As mentioned in chapter 1, this system 
automates much of the time-consuming and tedious task of exploring the behavior of 
dynamical systems. Typically, in order to characterize the behavior of a nonlinear cir¬ 
cuit, the investigator must select interesting parameter values with which to drive the 
circuit, measure the response, run a number of numerical computations (e.g. fast Fourier 
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transforms, autocorrelations and chirp z -transforms to help identify the behavior of the 
signal), and interpret the results. The results of these numeric computations are usually 
in the form of massive amounts of numerical data. Although the computer is very good 
at performing numeric computations, it is not very good at interpreting its qualitative 
content. Thus, it is largely up to the investigator to scan through the reams of data 
searching for relevant and important information. To have the analyzer program not 
only automatically set up and run the experiment, but produce a high-level qualitative 
description of the circuit’s behavior, clearly facilitates the task of the investigator. As 
mentioned previously, these high-level, qualitative descriptions are typically what an in¬ 
vestigator or experimental dynamidst is interested in. In fact, one of the aims of this 
project is to develop programs that can automatically generate descriptions of nonlinear 
dynamical systems similar to those found in published papere[AS87]. The content of the 
text from the high-level interpreter is based upon articles from the literature. 

6.1.2 Integration of instrument environment with a high-level pro¬ 
gramming language and environment 

This issue was touched upon briefly in chapter 2, but it is worth mentioning again here. 
The system provides a nice integration of the Scheme environment with the instrument 
environment. The advantage of such a combination is that it affords the investigator the 
flexibility of working in the Scheme environment, including such software packages as 
the signal-processing library, while also allowing access to a set of powerful measurement 
instruments. 

Another advantage of integrating the Scheme environment with the instrument envi¬ 
ronment is that it hides from the user the low-level details of how to communicate with 
and operate each instrument. This may not seem like a very large advantage if you only 
have one or two instruments to use. However, if you need to use many instruments, 1 each 
of which has a whole different sequence of buttons to press to achieve the desired effect, 
you begin to appredate a program which encapsulates this information. In our system, 
the communication module (refer to figure 2.3 of chapter 2) provides a nice interface be- 

1 The system for this thesis used six different instruments to explore the behavior of these nonlinear 
electrical circuits. 
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tween user and instruments. Thus, to perform a simple task like record-a-signal, the 
user need not open the owner’s manual and start reading about the two-million-and-one 
options on the waveform recorder. Instead, the analyzer program automatically knows 
how to select the appropriate parameter values for optimal operation of the instrument, 
and the communication module can send the appropriate, cryptic HP-IB command to 
the instrument for the desired result. This saves the user a great deal of time and grief, 
especially as the number of instruments used in the system increases. 

Finally, as mentioned in chapter 1, the integration of sophisticated measurement in¬ 
struments with powerful numerical packages and software programs, which interpret the 
results of physical measurements and numeric computations, suggests new possibilities 
for test and measurement. No longer are scientists limited to low-level, control-the-volt- 
meter type commands. They now have access to more powerful, more intelligent tools, 
capable of generating such high-level statements as, “The forced negative-resistance os¬ 
cillator undergoes a period doubling cascade as the driving frequency is increased from 
26 khz to 78 khz while the amplitude is fixed at 5.25 volts.” 

6.2 Future Work 

The Analyzer Program in combination with the rest of this system is a complex system. 
No doubt modifications and additions will be made to the system as it is tested on a 
greater number of nonlinear electrical circuits. However, I view this implementation as 
a successful first prototype. With regaxd to future work, I see several possibilities. 

One possible improvement to the system would be to increase the sophistication of 
the pattern recognizer in the high-level interpreter so it can recognize more sequences 
than it currently does. At present, the pattern recognizer is still at a fairly early stage of 
development and can only identify a few patterns. However, as we include more knowl¬ 
edge about common sequences in nonlinear dynamics, we can increase the number and 
complexity of patterns the high-level interpreter can recognize. Chapter 17 of Thompson 
and Stewart[TS86] provides a catalog of additional patterns we may wish to include in 
our knowledge base such as intermittency, the U sequence, alterna ting periodic-chaotic 
sequence, etc. 
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In addition to including more patterns in the knowledge base, we may wish to take ad¬ 
vantage of Feigenbaum’s universal number when searching for period-doubling cascades. 
M. J. Feigenbaum discovered a universal solution common to all systems undergoing 
period doubling. For a given system, if we denote by A n the value of the parameter at 
which its period doubles for the nth time, Feigenbaum’s discovery can be summarized 
by the equation 



where 6 n = 6 and 6 = 4.6692016... for large n[TS86, Fei83]. Thus, if we know the 
parameter values at which the period doubled the previous two times, we ran predict, 
approximately, the parameter value at which the period will double for the third time. 

Another improvement to the system might be to include a feedback loop. The ana¬ 
lyzer program, after performing a first pass analysis on a given region, should be able to 
automatically zoom in on regions of interest. For example, if the high-level interpreter 
sees al —► 2 —*• 4 —► 16 —*■ 32 —► ... —chaos sequence, the program should zoom 
into the region between the 4th-order subharmonic and the 16th-order subharmonic to 
search for the missing 8th-order subharmonic. In so doing, the program behaves in a 
more intelligent manner because it has some idea of what it is looking for; and if the 
program doesn’t find it, it looks more closely for the missing link. 

In addition to augmenting the knowledge base and including a feedback loop, the 
use of the autocorrelation technique to help distinguish between quasiperiodicity and 
chaos needs to be explored. As mentioned in chapter 4, the use of this technique looks 
promising. However, the details of implementation have yet to be worked out. 

For nonlinear systems, different initial conditions can lead to different steady state 
behaviors. Thus, to obtain truly reproducible parameter-space graphs, we need to ini¬ 
tialize the state variables of the nonlinear circuit to the same state each time. The best 
way to implement this is by adding additional hardware to the circuit. In the case of 
Ueda’s forced negative-resistance osdallator, for example, we mi gh t want add circuitry 
that automatically shorts out the capacitor and inductor prior to each test run. 

A final suggestion about improvements to the system would be to have the resolver 
module use knowledge about commonly recognized patterns, such as the period doubling 
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